A variable heart rate multi-compartmental coupled model of the cardiovascular and respiratory systems

Current mathematical models of the cardiovascular system that are based on systems of ordinary differential equations are limited in their ability to mimic important features of measured patient data, such as variable heart rates (HR). Such limitations present a significant obstacle in the use of such models for clinical decision-making, as it is the variations in vital signs such as HR and systolic and diastolic blood pressure that are monitored and recorded in typical critical care bedside monitoring systems. In this paper, novel extensions to well-established multi-compartmental models of the cardiovascular and respiratory systems are proposed that permit the simulation of variable HR. Furthermore, a correction to current models is also proposed to stabilize the respiratory behaviour and enable realistic simulation of vital signs over the longer time scales required for clinical management. The results of the extended model developed here show better agreement with measured bio-signals, and these extensions provide an important first step towards estimating model parameters from patient data, using methods such as neural ordinary differential equations. The approach presented is generalizable to many other similar multi-compartmental models of the cardiovascular and respiratory systems.


Introduction
Computational models of biological systems can help further understanding of physiology, and drive the development of predictive tools for clinical decisionmaking.A wide variety of models of the human cardiovascular system (CVS) have been proposed, and their complexity is frequently classified in terms of the number of spatial dimensions that are simulated [1,2].At the simplest end of the scale lie zero-dimensional (lumped parameter) models, which have no spatial information and instead model the CVS as discrete compartments (also known as a pressure-volume or PV model).One-dimensional models permit estimation of how flow and pressure waves propagate through a network of blood vessels [3][4][5].Two-dimensional, or more commonly, three-dimensional models use finite-element analysis (FEA) and/or computational fluid dynamics (CFD) to obtain highly detailed simulations of fluid flows [6,7]; these often still rely on lower dimensional models for the boundary conditions.
This work is concerned only with zero-dimensional time-dependent ordinary differential equation (ODE) models.Early examples of these models are the mono-compartment Windkessel [8,9] and Westkessel [10,11] models, employing electrical circuit analogies.By comparison, multi-compartment models split the CVS into discrete segments, enhancing the level of detail; using more compartments, however, requires more parameters to be estimated.One of the more extreme examples of this is the so-called Guyton model [12], which describes most of the main blood vessels as well as the regulation of autonomous and hormone systems.More commonly, simpler representations are used where the complexity can be increased in an area of interest [2,13,14].These models are typically only analysed once they reach a stable orbit, where the system returns to the same state after each heartbeat.As a result, such models are typically only simulated for as long as required for the initial transient behaviour (due to imprecise or unknown initial conditions) to settle; beyond this, no extra information can be obtained by running for longer.
The activity of the CVS and the respiratory system are strongly coupled, in part due to the location of the heart and adjacent arteries and veins inside the thoracic cavity, which experiences varying pressure due to activity of the respiratory muscles.Including these effects leads to more detailed models capable of investigating specific situations, such as the Valsalva manoeuvre or artificial ventilation, over longer time periods [2,[15][16][17][18][19][20][21].
These models may include some sort of regulatory system model for the control of heart rate (HR) and/or respiratory rate, leading to greatly increased model complexity [19,[21][22][23][24][25].Although this is useful for building understanding of specific conditions, implementing in silico all possible pathways by which respiration and circulation can be controlled in vivo is infeasible; for example, emotional responses can significantly affect HR, yet a model aiming to quantify or predict this in a clinical setting would be very challenging to design due to the plethora of underpinning mechanisms (and their quantification) that would be required.Some other contemporary CVS models, for example CircAdapt [26][27][28], permit open-loop control of HR, although it can only be changed in discrete steps and requires the ODE integration to be restarted at each change point (typically between predetermined 'rest' and 'exercise' settings).
This paper focuses on Smith's multi-compartment model of the human CVS with ventricular interaction [14,29], along with Jallon's extension to this model to simulate interaction between the CVS and the respiratory system [18].These models are designed to minimize complexity, and do not contain any mechanism for dynamic variation of HR.All models were implemented using JAX [30,31], with a view to permitting autodifferentiation of the ODE solutions.Future implications of this work, for example grey-box modelling, are discussed in §5.
Two novel modifications to these models are presented and explored in this paper.Firstly, a correction is derived in §2.3 to stabilize the behaviour of Jallon's respiratory system model over longer simulation times.Secondly, an extension to the cardiac driver function in the CVS model is developed in §2.4 to allow simulation incorporating a variable HR as an arbitrary function of time.
After specifying realistic model parameters and initial conditions in §2.5, results from simulations of the various models (including those previously published as well as those incorporating the extensions developed here) are presented in §3.The results demonstrate stable long-term cyclical behaviour of the interactive system even under variable HR conditions.A discussion follows in §4.
The approach to implementing variable HR developed herein is in theory generalizable to many other ODE models; the model given here is just one example of this.In §5, the paper concludes by discussing this, and also some further opportunities for applications of such models.

Methods
A review of Smith's model of the human CVS [14] is included in §2.1, along with an extension by Jallon [18] to simulate interaction between the heart and lungs in §2.2.Subsequently, §2.3 presents how Jallon's model can be stabilized to allow longterm simulation, and a method for implementing variable HR is developed in §2.4.Finally, §2.5 contains a complete listing of model parameters and initial conditions.
Python (JAX) implementations of the models described below are available on GitHub at slishak/cvsx.

Cardiovascular system model
The CVS has been minimally modelled as six compartments connected in a closed loop by valves, resistances and optionally inductances [14,22,32].A schematic is given in figure 1.This model has been widely applied in human and animal studies [1].It takes into account interaction between the ventricles through the ventricular septum and the pericardium.There are two models of valve dynamics that can be used: a simple noninertial check valve model which allows flow only during a negative (decreasing) pressure gradient, and a more realistic inertial model described as 'open on pressure, close on flow' which allows flow through the valve to continue due to the fluid inertia even if the pressure gradient becomes positive (increasing).royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230339 The model is described in full detail below, as there are a few cosmetic changes for readability compared with the referenced works (conditional statements are used instead of equivalent Heaviside functions).Each compartment's volume is governed by an ODE.

Single ventricle dynamics
The function of a ventricle in isolation is described by a PV diagram, which shows how the PV relationship changes through the cardiac cycle (figure 2).The two main characteristics are the end systolic pressure-volume relationship (ESPVR) and the end diastolic pressure-volume relationship (EDPVR).The ESPVR represents the maximum elastance of the ventricle at the end of systole (after contraction and ejection) and the EDPVR represents the minimum elastance during diastole (when the ventricle is relaxed and filling).Nomenclature of all variables and parameters in the equations below is given in tables 1 and 2, and all dimensional quantities are expressed in seconds, litres and mmHg except for HR which is expressed in min −1 by convention.
In this model, the equations governing the ESPVR and EDPVR, respectively, are given by where the subscript m specifies the virtual volume in accordance with the interactions between the ventricles and septum, as described below in §2.1.2[14].The full PV relationship applying to a ventricle, or time-varying elastance [33], is then given by where eðtÞ ¼ A i e ÀBiððtmodð60=HRÞÞÀCiÞ 2 : ð2:4Þ The function e(t) (equation (2.4)) is known as the cardiac driver function, which varies between 0 (diastole) and 1 (systole) over each cardiac cycle and is continuous and periodic with period 60/HR s; the multiplication by 60 is because HR is conventionally expressed in beats per minute, but t has units of seconds.A wide variety of such functions have been proposed, with complexities ranging from simple sums of Gaussian terms [34,35] to data based on invasive measurements of the CVS [36].The version presented in equation (2.4) is a clarified version of that from Smith's original model [14], with explicit modulo wrapping of t [1].Furthermore, note that when N = 1 (as commonly used), A 1 = 1 to satisfy 0 ≤ e(t) ≤ 1 and C 1 = 30/HR s (half a cardiac period) to satisfy continuity, leaving B 1 and HR as the only free parameters.A visualization of this choice of e(t) is shown in the bottom left of figure 3.An extension to the above to allow for variable HR is defined in §2.4.
ventricle pressure ventricle volume V 0 = volume at zero pressrure Example pressure-volume diagram for ventricular and septal walls (V lvf , V rvf , V spt ) over the cardiac cycle.© 2007 Elsevier Science & Technology Journals, from [22]; permission for use conveyed through Copyright Clearance Center, Inc.
Table 1.Pressure-volume relationships for the non-inertial cardiovascular system [22], with modifications after Jallon [18] in brackets (which are not used in any of the simulations below, but retained for reference).E es is the contractility, V d is the compartment dead-space, V 0 is the volume at zero pressure, and P 0 and λ are the gradient and curvature respectively for the EDPVR.The corresponding pressures are then calculated with equation (2.3) for m = lvf, rvf and spt, respectively.Because the pericardium sits inside the thorax, the thoracic ( pleural) pressure P pl is also applied (equation (2.8)), which is a constant in this model (but later becomes time-varying in §2.2 when covering Jallon's model).Similarly, the ventricles sit inside the pericardium so the pericardium pressure P peri (t) is applied to equations (2.9) and (2.10).The septal pressure is just the difference between the left and right ventricle pressure (equation (2.11)).Hence, To complete the system of equations, the septum volume V spt (t) must be found such that equations (2.3), (2.5), (2.6) and (2.11) are all satisfied.This does not have an analytic solution, but can be solved using a nonlinear root-finder applied to the resulting residual function, f ðV spt Þ ¼ P lvf ðV lv ðtÞ À V spt , tÞ À P rvf ðV rv ðtÞ þ V spt , tÞ À P spt ðV spt , tÞ: ð2:12Þ The gradient f 0 (V spt ) is also simple to compute by autodifferentiation or hand-calculation, which can be used for more efficient optimization (for example, using the Newton-Raphson algorithm).The uniqueness of the solution of f(V spt ) = 0 can be verified by observing from equation (2.3) that dP m /dV m > 0 for any set of physical ( positive) parameter values, and therefore that f 0 (V spt ) is negative for any V spt value (and so f(V spt ) cannot cross zero more than once).Note that linear approximations to the system above have been proposed [18,37], but they have generally been found here to be unnecessary if root finder is initialized with the solved V spt from the last accepted ODE step, as the Newton-Raphson algorithm then converges to an acceptable tolerance within one or two iterations.However, when trying to backpropagate gradients through the ODE solution [38], specialized solvers [31,39] may be required.Furthermore, this backpropagation through a nonlinear solver can represent an additional slowdown.With this in mind, Jallon's EDPVR linearization [18], given by P ed,lin,spt ðVÞ ¼ P 0,spt l spt ðV À V 0,spt Þ, ð2:13Þ replaces equation (2.2) only in the solution for V spt , yielding the full form of which has not been explicitly given before.In equation (2.14), the (t) argument has been hidden from V spt (t), V lv (t) and V rv (t) to reduce clutter.This optional linearization can result in substantially different behaviour, shown later in this paper.

Valve dynamics
Neglecting inertial effects, the equation for the instantaneous flow rate Q v (t) induced within an internal vasculature v of resistance R v from an upstream pressure P up,v (t) to a downstream pressure P down,v (t) is given by where L v is a constant inductance [40].An 'open on pressure, close on flow' valve law is used.Note that as long as Q v (0) is Table 2. Other parameters for the non-inertial cardiovascular system [22], with modifications after Jallon [18]  positive, in theory the valve law as stated should prevent the flow rate ever becoming negative (as dQ v /dt can never be negative when Q v (t) = 0), but in practice, the tolerances in ODE solvers mean that there is typically some small constant negative flow rate; this is quoted as typically being between −1 × 10 −4 and −1 × 10 −6 [37] but depends on the choice of ODE solver used, the units of the states, and the error (step size) control.The model as presented here is equivalent to the full Heaviside formulation in [37].Regardless of whether inertia is being considered, the rate of change volume of a single chamber c is then described by equation (2.19), assuming incompressible fluid: the rate of change of volume (mass) in the chamber is equal to the net inflow of volume (mass).To account for the tiny negative flow rate that is possible when considering inertia (equation (2.18)), the ramp function in equation (2.17) is used to prevent the error propagating around the rest of the model (although it is not strictly necessary with the non-inertial model, as it is already applied in equation (2.15)), for c = rv, pa, pu, lv, ao and vc, and the flow rates in/out correspond to those shown in figure 1 for compartment c.The valve conditions shown in equations (2.17) and (2.18) are also suitable for use with the event-handling capabilities of some ODE solvers, by only switching from one branch to the other when an exact zero crossing of the condition is found.As this V spt (ml) e (t) Smith non-inertial and inertial CVS model outputs with parameters from tables 3 and 4, providing a like-for-like comparison on the effects of including inertia in the model.
is not currently supported by the chosen ODE solver [31], this might be a future improvement, although the use of an adaptive step ODE solver already achieves a low error around this discontinuity [37,41].

Full closed-loop model
The final model is formed by setting the flow rate out of each compartment in figure 1 equal to the flow rate into the next compartment.The implementation used to generate the results below is almost the same as [1,37] but with the pleural pressure P pl applied to the pulmonary vein and artery pressures P pa and P pu [18].The total blood volume, is conserved throughout the simulation, as it is a closed fluid flow system.The initial individual compartment volumes are less important, as they will not affect the asymptotic behaviour of the system as long as they are feasible.Stable initial compartment volumes can be found manually by running the model from a reasonable starting point (e.g.tables 6 or 7) and recording a set of state values from the orbit that the model converges to.Some other possible approaches are discussed in [29].Initial flow rates for the inertial model can be found using equation (2.15) (i.e.assuming zero inertia).

Combining respiratory and cardiovascular models
The CVS model presented above stabilizes to periodic behaviour, with every heartbeat exhibiting the same dynamics.However, in reality, the pleural pressure P pl is not constant but varies with the respiratory rate, which (among many other factors) causes beat-to-beat variation in the cardiovascular dynamics.In Jallon's heart-lung model [18], the respiratory rate is controlled by the central respiratory pattern generator, which is based on the Liènard system [42]   The full respiratory/cardiovascular model can be simulated as a system of ODEs, constructing a right-hand side function defining the state derivatives using the following procedure (after setting appropriate initial conditions and parameter values, discussed in §2.5): (i) receive time t and states (CVS volumes and flow rates, x(t), y(t), P mus (t) and V alv (t)), (ii) call the respiratory system model (equations ( 2

Correcting and stabilizing the respiratory model
Equation (2.24) in §2.2 has been found to cause a gradual drift of P mus over time, resulting in long-term instability of the solution and unbounded increase or decrease in alveolar volume V alv .This is visible in Jallon's original paper (figure 5) and is also possibly the cause of issues from other referenced uses of this model [45].A later paper from the same author introduces a slightly modified central respiratory pattern generator model, but this is understood to be a simplification in order to simulate a specific phenomenon rather than a general model improvement (and parameter unit inconsistencies further confuse the matter) [46].Although equation (2.21) contains a term that supposedly prevents over-inflation of the lungs through the Hering-Breuer reflex, the justification for this implementation is unclear (an increase in V alv acts to decrease the third derivative of P mus ).Although these reflexes were historically thought to play a significant role in the regulation of ventilation, they are now understood to be largely inactive in adult humans (except while exercising), although there is some evidence that they may be more important in newborn babies [47].Regardless, as the respiratory pattern generator is only required to synthesize realistic behaviour for P mus , there is little need for the complication of trying to mimic regulatory systems in the model (and furthermore, no clear reason for implementing just one element of a complicated control system).
Below, a modification is proposed to add an extra term to essentially perform integral feedback control (a type of proportional-integral-derivative or PID control [48]) on P mus , introducing a new parameter β in equation (2.29) (replacing equation (2.24)) which stabilizes the solution for small positive β (for example, β = 0.1 s −1 ).As this acts to prevent P mus and V alv steadily increasing or decreasing over time, the attempt to implement the Hering-Breuer reflex in equation ( 2  royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230339

Model parameters and initial conditions
Full parametrizations for all biophysical models are included in tables 1 to 5. As detailed below, the parametrization is taken from previously published works, except where they are believed to be erroneous, as the contribution of this work is methodological in nature and does not depend on any specific novel parametrization.All parameters are converted into consistent units of pressure (mmHg) and volume (litres).
Complete and accurate parameter listings for the inertial cardiovascular model are rare in the existing published literature due to the challenges of in situ measurement and variations across animal models, species etc.One useful source is [1], but with a few necessary modifications: the given value for E ao is 0, which is assumed to be a mistake as it results in non-physical behaviour, and R pul and R sys are missing.The values in tables 3 and 4 came from a later PhD thesis [49].Furthermore, the referenced value of the offset parameter C i in the cardiac driver function (equation (2.4)) appears to be incorrect, as setting it to anything other than half the cardiac cycle duration results in a discontinuity at the end of the cycle.Finally, V tot is described as 5.5 l in the referenced source, but due to the low total deadspace volume in table 3, this results in non-physical behaviour, therefore a stressed blood volume of 1.5 l is assumed (and the model ignores the 4 l of dead space).Another solution would have been to restore the dead space in the PV compartments, which sum to 4 l (excluding the virtual septum volume).
Table 5 gives Jallon's parameters for the respiratory model.In tables 1 and 2, there are some alternative parameter values provided by Jallon [18] in brackets which reduce the stiffness of the septum; these are described as being necessary to achieve reasonable physiological values for all variables.However, no physical justification was made for this change, and they cause the ventricle PV loops (figure 2) to lose their shape completely.As a result, these parameters are not used in the simulations.
All initial state values are given in tables 6, 7 and 8.The initial compartment volumes must sum to V tot (and this is the only place the parameter is used).Apart from this, the asymptotic behaviour of the model does not depend on the individual initial volumes, as long as they are sensible.The closer the initialization is to the asymptotic behaviour, the shorter the initial transient period will be.Alternatively, initialization could be performed by an optimization procedure [29].
The models were simulated using Diffrax [31] with the Tsit5 solver, frequently recommended as a good generalpurpose ODE solver [31,50,51].Adaptive step sizing was used, with the default proportional error control; this is important Table 3. Pressure-volume relationships for the inertial cardiovascular system [1].Parameters marked with (*) differ from the source, for reasons explained in the main text.1088.9 mmHg s l −1 (*) total blood volume (V tot ) 1.5 l (*) Table 5. Parameters for respiratory model [18].Note that R ua , R ca and V th , 0 are given with incorrect units in the original paper, and HB, λ and μ are given with no units.The parameter β is introduced in §2.3.Bear in mind that the Liènard system parameters are distinct from the similarly named A i , B i in the cardiac driver function.
The absolute tolerance was set to 10 −7 as the unwanted negative flow rates described in §2.1.3were found to respond linearly to this with a gain of roughly 100, giving typical values of around −10 −5 1 s −1 .The relative tolerance was set to 10 −4 , as with the default setting of 10 −3 , the interpolated dense solution for the non-inertial model was noticeably inaccurate around the closing of the pulmonary and aortic valves, with oscillating flow rates between ODE steps.

Results
Results of simulating the various models presented in §2 are shown below.The dense ODE solution is interpolated onto a regular grid, so that a maximum step size does not need to be set just for visualization purposes.States were expressed in units of litres (but converted to millilitres for plotting), and pressures were calculated in mmHg.Figure 3 shows two seconds of simulation for the inertial and non-inertial Smith CVS models with ventricular interaction, using parameters from tables 3 and 4 in both cases.No previous published comparison of these two models with matching parameters has been found.In the upper right corner of the PV diagram, the inertial model is less rounded due to the sharper ejection profile, which is also visible in Q av and Q pv .Other than this, both models are very similar; in particular, the septum deviation volume V spt is almost unaffected.Figure 6 shows the behaviour of the Jallon heart-lung model presented in §2.2 [18].Following the original implementation, the linearization of V spt is used, although the suggested parameter modifications (such as reducing the stiffness of the septum) were discarded, as explained in §2.5.The behaviour matches expectations, with left ventricle stroke volume decreasing during inspiration and increasing during expiration, and vice versa for the right ventricle.However, a gradual drift in almost all outputs can be seen, due to a growth of P mus over time.
Figure 7 gives a clear comparison of the effects of Jallon's linearization of the ventricular interaction.Although the linearization has little influence on the numerical results during periods of static pleural pressure (or when the lungs are mostly exhaled), there is a noticeable difference during diastole immediately after inhalation.The full nonlinear model suggests a substantial septal deflection into the left ventricle in this situation, but the linearized model does not show this.
The results of the correction proposed in §2.3 that intends to stabilize Jallon's respiratory model are shown in figure 8, confirming that there is no longer any drift over time.The HB parameter has been set to zero as the reflex has been accounted for by setting β = 0.1.The linearized ventricular interaction model was retained for a direct comparison against the original model.
Figure 9 shows the results of the inertial Smith model with the variable HR extension proposed in §2.4.As the HR gradually increases from 60 up to 100, following HRðtÞ ¼ 80 þ 20 tanh ð0:3ðt À 20ÞÞ, ð3:1Þ the peak flow rates through the aortic and pulmonary valves increase, and systolic/diastolic arterial pressures all increase substantially.The full nonlinear V spt solver (using equation (2.12)) was used in this simulation.Table 9 shows the computation time achieved in JAX for each of the model configurations simulated.

Discussion
It is difficult to evaluate the accuracy of the results presented in §3 in comparison with the published versions, due to the lack of provided code and occasionally insufficient descriptions of the parameters.However, qualitative comparison with the available plots suggests that all models were correctly reproduced.Upon inhalation, right ventricle stroke volume increases and left ventricle decreases, as expected.
The V spt linearization defined by equation (2.14) had a much larger than expected effect on the V spt trace when combined with the Jallon model.During inspiration, it results in a substantial under-prediction of septal deflection compared with the full nonlinear model.Recall that the linearization is motivated by computational efficiency rather than accuracy, so regardless of which is more physically accurate, the intention of the model has been substantially affected.However, the criteria that the V spt variation should be around 4% of the ventricle, quoted in Smith's original work [14,34], is no longer met in the Jallon model with nonlinear ventricular interaction, and the right ventricle expands notably.If the septal elasticity were further reduced, as suggested by Jallon, this effect would be even more pronounced.The computational advantages of the linearized ventricular interaction were found here to be limited.In an earlier PyTorch implementation, the nonlinear Newton-Raphson solver already converged within one or two iterations when initialized to the value of the last accepted ODE step (so the improvement was negligible).In the final JAX implementation, varying the initial guess was not found to be possible.Despite this, even though the linearization resulted in a speed-up of over three times, the model was already capable of simulating one second of activity in under 1 ms.However, the linearization may be useful in a parameter estimation setting where the simulation might need to be iterated, and the backwards pass may also be sped up with the linearization.
The modification proposed in §2.3 to stabilize the respiratory model is considered to be successful, and is necessary for long-term simulations of the Jallon model.The modification does not affect the physical relevance of the model, as it only slightly adjusts the use of the Liènard system as a signal generator.
The variable HR model proposed in §2.4 also works as intended, allowing for ICU HR data to be fed directly into the CVS models.This change would be useful when performing parameter estimation on long periods of patient data, as it provides an alternative to the beat-to-beat approach [52].On the other hand, the physiological relevance of the model results remain to be verified.Figure 9 suggests that when the HR increases, the diastolic aortic pressure should noticeably increase, with the systolic pressure also increasing but at a shallower rate.However, studies on patients fitted with pacemakers suggest that in humans, although the diastolic aortic pressure does indeed increase when the HR is artificially raised (as do peripheral systolic and diastolic arterial pressures), the systolic aortic pressure does not increase [53,54].This effect is understood to be due to wave reflection and change in systemic arterial stiffness, mechanisms that are not simulated within this model framework.Such neglected aspects go some way to explaining the difference in behaviour between the model and experimental results.royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230339 There are a wide range of mechanisms for the autonomous regulation of blood pressure [55] that are not considered by this model.Furthermore, an assumption is (most likely incorrectly) made that the length of systole as a proportion of the cardiac period does not change with HR, even when the driving pulse (equation (2.33)) is squashed or stretched in time depending on the rate of change of s w (t).Finally, the shape of the cardiac driver function used within this model is extremely simplistic, being defined by a single parameter alone ( B1 ).It would in theory be possible to redefine this parameter to be a function of HR; this approach could also be extended to more complicated driver functions [34][35][36]56].
Although it is not demonstrated in this work, the use of an ODE solver that supports autodifferentiation may permit the parameters to be learned directly by minimizing a loss versus clinically measured patient-specific blood pressure and HR data.Furthermore, it could lead to a grey-box model that combines human-designed models of the CVS with regulated variables such as HR, elastances and resistances implemented as unknown functions (of time and/or state variables) to be learned, for example neural ODEs [38].Such models would be valuable as they could estimate information about a patient that would be infeasible or overly intrusive to measure directly.In particular, this approach of combining a simple model with learned functions would address the issue of increasingly complex models having many unknown physiological parameters which cannot be uniquely identified from commonly available patient data.In conjunction with real-time data available in an ICU, this could provide vital extra information when making important decisions on treatments.royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230339

Conclusion
Two novel modifications to the CVS models of Smith [14] and Jallon [18] were proposed: a correction to stabilize Jallon's model ( §2.3) and an extension to allow all models to simulate a variable HR ( §2.4).Both modifications work as intended and significantly extend the capabilities of Smith's original model.The stabilization of the Jallon model is crucial to enable realistic simulation of, and parameter estimation from, physiological behaviour over clinically relevant time scales, and supersedes the previous attempt to incorporate the Hering-Breur reflex into the original model.
The novel approach to simulating variable HR, although applied here on a very simple model and basic cardiac driver function, has also been successful.However, the model suggests an increase in central systolic arterial pressure with increasing HR, which is not in agreement with available experimental data.It is postulated that this discrepancy may be due to the model neglecting the effects of pressure wave reflection and changes in systemic arterial stiffness with HR.On the other hand, the extension is generalizable to many other ODE models of the CVS, so a similar modification to a more complicated original model might yield more accurate results.For example, the CircAdapt [26] model might also be extended in this manner.V spt (ml) ventricle pressure (mmHg) royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230339 With the HR defined by an arbitrary function of time, the model is no longer limited to periodic behaviour.In conjunction with the support for autodifferentiation due to the implementation in JAX, this leads to further opportunities for applying the CVS model in conjunction with machinelearning techniques to measured bio-signals from patients, especially when considering the high computational cost of parameter estimation with models which do not support differentiation and have higher parameter counts [57].
Ethics.This work did not require ethical approval from a human sub- ject or animal welfare committee.

Figure 1 .
Figure 1.Schematic of the full cardiovascular system model.The inductances are only present when modelling inertial flows.The diodes represent check valves.Nomenclature is defined in tables 1 and 2.

Figure 4 .
Figure 4. Diagram explaining how the left and right ventricles are split into three free wall volumes inside the pericardium.© 2005 Elsevier Science & Technology Journals, from [32]; permission for use conveyed through Copyright Clearance Center, Inc.

Figure 5 .
Figure 5. Plots of P mus and V alv from the Jallon heart-lung model, which show undesired non-stationary behaviour with P mus steadily increasing and V alv decreasing.© 2009 The Royal Society (UK), from [18]; permission for use conveyed through Copyright Clearance Center, Inc. (with overlaid red lines to accentuate the drift).

Figure 6 .Figure 7 .
Figure 6.Selected outputs from the Jallon heart-lungs model, with parameters from tables 1, 2 and 5.The original behaviour (β = 0) is plotted in blue, with the stabilized model (β = 0.1) overlaid in red.Note the gradual drift in behaviour.

Figure 9 .
Figure 9. Smith inertial cardiovascular system model outputs with parameters from tables 3 and 4. Variable heart rate simulated with model from §2.4 and equation (3.1).

Table 4 .
[1]er parameters for the inertial cardiovascular system[1].Parameters marked with (*) differ from the source, for reasons explained in the main text.The cardiac driver and thoracic pressure parameters are as table2.

Table 6 .
Initial volumes for non-inertial cardiovascular system model parametrized by tables 1 and 2.

Table 7 .
Initial volumes for inertial cardiovascular system model parametrized by tables 3 and 4. Initial flow rates can be calculated using equation(2.15).Note that the unstressed blood volume is not considered in this parametrization.

Table 8 .
Initial volume ratios for respiratory model parametrized by table 5.The Liènard system states are both dimensionless.

Table 9 .
Computation time in JAX for 60 s of simulation (not including compilation time), averaged over 10 runs.Run on an AMD Ryzen 7 5800H CPU in Windows Subsystem for Linux.Linear/nonlinear refers to ventricular interaction; the Jallon model is always non-inertial.The compiled-in max number of steps in Diffrax was set to 164..org/journal/rsifJ. R. Soc.Interface 20: 20230339 royalsocietypublishing